/*
 copyright 2020, Chen Wei <weichen302@gmail.com>
 version 0.0.3
Implement astronomical algorithms for finding solar terms and moon phases.

Truncated IAU2000B from USNO NOVAS c source is used for nutation.

The higher accuracy light abberation algorithm from A & A. p156,

*/

#include <stdio.h>
#include <math.h>
#include "astro.h"

/*
 77 terms IAU2000B Luni-Solar nutation table from USNO NOVAS c source

 Luni-Solar argument multipliers,  coefficients, unit 1e-7 arcsec
 L  L'  F  D  Om longitude (sin, t*sin, cos), obliquity (cos, t*cos, sin)
*/
static double IAU2000BNutationTable[77][11] = {
    {  0,  0,  0,  0, 1, -172064161, -174666,  33386, 92052331,  9086, 15377 },
    {  0,  0,  2, -2, 2,  -13170906,   -1675, -13696,  5730336, -3015, -4587 },
    {  0,  0,  2,  0, 2,   -2276413,    -234,   2796,   978459,  -485,  1374 },
    {  0,  0,  0,  0, 2,    2074554,     207,   -698,  -897492,   470,  -291 },
    {  0,  1,  0,  0, 0,    1475877,   -3633,  11817,    73871,  -184, -1924 },
    {  0,  1,  2, -2, 2,    -516821,    1226,   -524,   224386,  -677,  -174 },
    {  1,  0,  0,  0, 0,     711159,      73,   -872,    -6750,     0,   358 },
    {  0,  0,  2,  0, 1,    -387298,    -367,    380,   200728,    18,   318 },
    {  1,  0,  2,  0, 2,    -301461,     -36,    816,   129025,   -63,   367 },
    {  0, -1,  2, -2, 2,     215829,    -494,    111,   -95929,   299,   132 },
    {  0,  0,  2, -2, 1,     128227,     137,    181,   -68982,    -9,    39 },
    { -1,  0,  2,  0, 2,     123457,      11,     19,   -53311,    32,    -4 },
    { -1,  0,  0,  2, 0,     156994,      10,   -168,    -1235,     0,    82 },
    {  1,  0,  0,  0, 1,      63110,      63,     27,   -33228,     0,    -9 },
    { -1,  0,  0,  0, 1,     -57976,     -63,   -189,    31429,     0,   -75 },
    { -1,  0,  2,  2, 2,     -59641,     -11,    149,    25543,   -11,    66 },
    {  1,  0,  2,  0, 1,     -51613,     -42,    129,    26366,     0,    78 },
    { -2,  0,  2,  0, 1,      45893,      50,     31,   -24236,   -10,    20 },
    {  0,  0,  0,  2, 0,      63384,      11,   -150,    -1220,     0,    29 },
    {  0,  0,  2,  2, 2,     -38571,      -1,    158,    16452,   -11,    68 },
    {  0, -2,  2, -2, 2,      32481,       0,      0,   -13870,     0,     0 },
    { -2,  0,  0,  2, 0,     -47722,       0,    -18,      477,     0,   -25 },
    {  2,  0,  2,  0, 2,     -31046,      -1,    131,    13238,   -11,    59 },
    {  1,  0,  2, -2, 2,      28593,       0,     -1,   -12338,    10,    -3 },
    { -1,  0,  2,  0, 1,      20441,      21,     10,   -10758,     0,    -3 },
    {  2,  0,  0,  0, 0,      29243,       0,    -74,     -609,     0,    13 },
    {  0,  0,  2,  0, 0,      25887,       0,    -66,     -550,     0,    11 },
    {  0,  1,  0,  0, 1,     -14053,     -25,     79,     8551,    -2,   -45 },
    { -1,  0,  0,  2, 1,      15164,      10,     11,    -8001,     0,    -1 },
    {  0,  2,  2, -2, 2,     -15794,      72,    -16,     6850,   -42,    -5 },
    {  0,  0, -2,  2, 0,      21783,       0,     13,     -167,     0,    13 },
    {  1,  0,  0, -2, 1,     -12873,     -10,    -37,     6953,     0,   -14 },
    {  0, -1,  0,  0, 1,     -12654,      11,     63,     6415,     0,    26 },
    { -1,  0,  2,  2, 1,     -10204,       0,     25,     5222,     0,    15 },
    {  0,  2,  0,  0, 0,      16707,     -85,    -10,      168,    -1,    10 },
    {  1,  0,  2,  2, 2,      -7691,       0,     44,     3268,     0,    19 },
    { -2,  0,  2,  0, 0,     -11024,       0,    -14,      104,     0,     2 },
    {  0,  1,  2,  0, 2,       7566,     -21,    -11,    -3250,     0,    -5 },
    {  0,  0,  2,  2, 1,      -6637,     -11,     25,     3353,     0,    14 },
    {  0, -1,  2,  0, 2,      -7141,      21,      8,     3070,     0,     4 },
    {  0,  0,  0,  2, 1,      -6302,     -11,      2,     3272,     0,     4 },
    {  1,  0,  2, -2, 1,       5800,      10,      2,    -3045,     0,    -1 },
    {  2,  0,  2, -2, 2,       6443,       0,     -7,    -2768,     0,    -4 },
    { -2,  0,  0,  2, 1,      -5774,     -11,    -15,     3041,     0,    -5 },
    {  2,  0,  2,  0, 1,      -5350,       0,     21,     2695,     0,    12 },
    {  0, -1,  2, -2, 1,      -4752,     -11,     -3,     2719,     0,    -3 },
    {  0,  0,  0, -2, 1,      -4940,     -11,    -21,     2720,     0,    -9 },
    { -1, -1,  0,  2, 0,       7350,       0,     -8,      -51,     0,     4 },
    {  2,  0,  0, -2, 1,       4065,       0,      6,    -2206,     0,     1 },
    {  1,  0,  0,  2, 0,       6579,       0,    -24,     -199,     0,     2 },
    {  0,  1,  2, -2, 1,       3579,       0,      5,    -1900,     0,     1 },
    {  1, -1,  0,  0, 0,       4725,       0,     -6,      -41,     0,     3 },
    { -2,  0,  2,  0, 2,      -3075,       0,     -2,     1313,     0,    -1 },
    {  3,  0,  2,  0, 2,      -2904,       0,     15,     1233,     0,     7 },
    {  0, -1,  0,  2, 0,       4348,       0,    -10,      -81,     0,     2 },
    {  1, -1,  2,  0, 2,      -2878,       0,      8,     1232,     0,     4 },
    {  0,  0,  0,  1, 0,      -4230,       0,      5,      -20,     0,    -2 },
    { -1, -1,  2,  2, 2,      -2819,       0,      7,     1207,     0,     3 },
    { -1,  0,  2,  0, 0,      -4056,       0,      5,       40,     0,    -2 },
    {  0, -1,  2,  2, 2,      -2647,       0,     11,     1129,     0,     5 },
    { -2,  0,  0,  0, 1,      -2294,       0,    -10,     1266,     0,    -4 },
    {  1,  1,  2,  0, 2,       2481,       0,     -7,    -1062,     0,    -3 },
    {  2,  0,  0,  0, 1,       2179,       0,     -2,    -1129,     0,    -2 },
    { -1,  1,  0,  1, 0,       3276,       0,      1,       -9,     0,     0 },
    {  1,  1,  0,  0, 0,      -3389,       0,      5,       35,     0,    -2 },
    {  1,  0,  2,  0, 0,       3339,       0,    -13,     -107,     0,     1 },
    { -1,  0,  2, -2, 1,      -1987,       0,     -6,     1073,     0,    -2 },
    {  1,  0,  0,  0, 2,      -1981,       0,      0,      854,     0,     0 },
    { -1,  0,  0,  1, 0,       4026,       0,   -353,     -553,     0,  -139 },
    {  0,  0,  2,  1, 2,       1660,       0,     -5,     -710,     0,    -2 },
    { -1,  0,  2,  4, 2,      -1521,       0,      9,      647,     0,     4 },
    { -1,  1,  0,  1, 1,       1314,       0,      0,     -700,     0,     0 },
    {  0, -2,  2, -2, 1,      -1283,       0,      0,      672,     0,     0 },
    {  1,  0,  2,  2, 1,      -1331,       0,      8,      663,     0,     4 },
    { -2,  0,  2,  2, 2,       1383,       0,     -2,     -594,     0,    -2 },
    { -1,  0,  0,  0, 2,       1405,       0,      4,     -610,     0,     2 },
    {  1,  1,  2, -2, 2,       1290,       0,      0,     -556,     0,     0 }
};

/*
    Calculate nutation angles using the IAU 2000B model.

    The table is from NOVAS(USNO) c source flle.

    "...the 77-term, IAU-approved truncated nutation series, IAU 2000B, which
    is accurate to about 0.001 arcsecond in the interval 1995-2050"

    Arg:
        jd as jd
    Return:
        nutation of longitude in radians
 */
double nutation(double jd) {
    double t, L, Lp, F, D, Om;
    t = (jd - J2000) / 36525.0;

    /* Mean anomaly of the Moon, in arcsec */
    L = 485868.249036 + t * 1717915923.2178;

    /* Mean anomaly of the Sun. */
    Lp = 1287104.79305 + t * 129596581.0481;

    /* Mean argument of the latitude of the Moon. */
    F = 335779.526232 + t * 1739527262.8478;

    /* Mean elongation of the Moon from the Sun. */
    D = 1072260.70369 + t * 1602961601.2090;

    /* Mean longitude of the ascending node of the Moon. */
    Om = 450160.398036 - t * 6962890.5431;

    double arg, lon, dpplan;
    lon = 0;
    int i;
    for (i = 0; i < 77; i++) {
        arg = (   IAU2000BNutationTable[i][0] * L
                + IAU2000BNutationTable[i][1] * Lp
                + IAU2000BNutationTable[i][2] * F
                + IAU2000BNutationTable[i][3] * D
                + IAU2000BNutationTable[i][4] * Om) * ASEC2RAD;

        lon +=    (  IAU2000BNutationTable[i][5]
                   + IAU2000BNutationTable[i][6] * t) * sin(arg)
                + IAU2000BNutationTable[i][7] * cos(arg);
    }

    /* unit of longitude is 1.0e-7 arcsec, convert it to arcsec */
    lon *= 1.0e-7;

    /* Constant account for the missing long-period planetary terms in the
       truncated nutation model, in arcsec */
    dpplan = -0.000135;
    lon += dpplan;

    lon *= ASEC2RAD;
    /* lon = fmod(lon, TWOPI); */
    return lon;
}

/* the higher accuracy light abberation algorithm from A & A. p156,
 * the error will be less than 0.001"
 */
double lightabbr_high(double jd)
{
    double tm, tm2, tm3, t, t2, t3;
    t = (jd - J2000) / 36525.0;
    t2 = t * t;
    t3 = t * t2;

    tm = t / 10.0;
    tm2 = tm * tm;
    tm3 = tm * tm2;

    // variation of the Sun's longitude
    double var_lon;
    var_lon = 3548.330
        + 118.568 * sin(1.527664004990 +   6283.075850600876 * tm)
        +   2.476 * sin(1.484508993906 +  12566.151699456424 * tm)
        +   1.376 * sin(0.486077687339 +  77713.771468687730 * tm)
        +   0.119 * sin(1.276490181677 +   7860.419392621536 * tm)
        +   0.114 * sin(5.885711004647 +   5753.384884566103 * tm)
        +   0.086 * sin(3.884055717388 +  11506.769769132206 * tm)
        +   0.078 * sin(2.841633387025 + 161000.685738008644 * tm)
        +   0.054 * sin(1.441333038870 +  18849.227550057301 * tm)
        +   0.052 * sin(2.993569534399 +   3930.209696310768 * tm)
        +   0.034 * sin(0.529208263814 +  71430.695618086858 * tm)
        +   0.033 * sin(2.091087703461 +   5884.926847413107 * tm)
        +   0.023 * sin(4.320419446313 +   5223.693920276658 * tm)
        +   0.023 * sin(5.674983441420 +   5507.553238331428 * tm)
        +   0.021 * sin(2.707426294193 +  11790.629088932305 * tm)
        +   7.311 * tm  * sin(5.819826570714 +   6283.075850600876 * tm)
        +   0.305 * tm  * sin(5.776715192860 +  12566.151699456424 * tm)
        +   0.010 * tm  * sin(5.733703298774 +  18849.227550057301 * tm)
        +   0.309 * tm2 * sin(4.214128894867 +   6283.075850600876 * tm)
        +   0.021 * tm2 * sin(3.578766215288 +  12566.151699456424 * tm)
        +   0.004 * tm3 * sin(5.198655163283 +  77713.771468687730 * tm)
        +   0.010 * tm3 * sin(2.700139544566 +   6283.075850600876 * tm);

    double M, e, C, v, R, res;
    // mean anomaly of the Sun
    M = (357.52910 + 35999.0503 * t - 0.0001559 * t2
                                    - 0.00000048 * t3) * DEG2RAD;
    // the eccentricity of Earth's orbit
    e = 0.016708617 - 0.000042037 * t - 0.0000001236 * t2;
    // Sun's equation of center
    C = (  (1.9146 - 0.004817 * t - 0.000014 * t2) * sin(M)
         + (0.019993 - 0.000101 * t) * sin(2 * M)
         + 0.00029 * sin(3 * M)) * DEG2RAD;
    // true anomaly
    v = M + C;
    // Sun's distance from the Earth, in AU
    R = (1.000001018 * (1 - e * e)) / (1 + e * cos(v));

    res = -0.005775518 * R * var_lon * ASEC2RAD;

    return res;
}

